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A brief summary of the formulation of QCD at finite chemical potental, /i, is presented. The failure of the 
quenched approximation to the problem is reviewed. 

Results are presented for dynamical simulations of the theory at strong and intermediate couplings. We find 
that the problems associated with the quenched theory persist: the onset of non-zero quark number does seem 
to occur at a chemical potential ~ However analysis of the Lee- Yang zeros of the grand canonical partition 
function in the complex fugacity plane, (e M,/T ), does show signals of critical behaviour in the expected region of 
chemical potential. 

Results are presented for a simulation at finite density of the Gross-Neveu model on a 16 3 lattice near to the 
chiral limit. Contrary to our simulations of QCD no pathologies were found when fj, passed through the value 
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1. Introduction 

Lattice QCD at finite baryon density holds the 
key to an understanding of the phase transition 
between quark-confined hadronic matter and the 
quark-gluon plasma. At zero temperature we ex- 
pect this transition to occur at [i = ^ corre- 
sponding to the lowest lying state with non-zero 
baryon number. 

Although the QCD phase transition at finite tem- 
perature has been successfully studied on the 
lattice with simulations predicting a first order 
transition at a well defined critical temperature, 
a non-perturbative lattice calculation of nuclear 
matter at finite density has been an outstanding 
problem in QCD thermodynamics since 1986 @,|. 
The fundamental difficulty in simulating QCD 
at finite density and investigating the transition 
quantitatively is that the effective action resulting 
from the Grassmann integration over the fermions 



is complex due to the introduction of the chemi- 
cal potential in the Dirac matrix. This complex 
nature of the QCD finite density action |?]J4|] pro- 
hibits the use of naive probabilistic methods in 
evaluating the functional integral. Thus the stan- 
dard Hybrid Monte Carlo simulation algorithms 
for lattice QCD with dynamical fermions are in- 
appropriate in this context since they require a 
positive definite probability measure. 

The study of lattice QCD at finite density is 
motivated by the need for a full equation of state 
for nuclear matter in the temperature-chemical 
potential plane which would give quantitative 
predictions for critical energies and densities in 
relativistic heavy ion collisions. Other applica- 
tions include equations of state for neutron stars 
and cosmological models of the early universe. 



2 



2. Quenched Simulations 

Serious problems were first reported in 
quenched simulations of finite density QCD in 
1986|Q] and the physical and mathematical rea- 
sons for this failure have been the focus of con- 
siderable debate ever since @||. 

In these early simulations the behaviour of the 
chiral condensate was studied at fixed quark mass 
for various different values of the chemical poten- 
tial, /i. The behaviour was initially as expected 
viz. the chiral condensate remained constant up 
until a certain value of mu and then tended to 
zero as the chemical potential was increased. 

We would expect that physical observables are 
/i independent up to some fi c which is related to 
the threshold for baryon production. The prob- 
lem encountered was the following: as the quark 
mass was decreased the \i at which the chiral con- 
densate began to change also decreased towards 
zero. 

In fact, the onset of chiral symmetry restora- 
tion appeared to occur at a chemical potential of 
half the pion mass, which would extrapolate to 
zero in the chiral limit. A study of the distribu- 
tion of the eigenvalues of the lattice Dirac opera- 
tor confirmed that the lowest mass state contain- 
ing a net number of fermions (i.e. a quark or a 
baryon but not a meson which does not see \x) 
became massless as the bare quark mass was re- 
duced to zero. The interpretation was that there 
existed either baryonic states which became mass- 
less in the chiral limit and had an energy equal 
to \ra^ or stable quark matter with low mass per 
baryon. 

The result that fj, c is proportional to the pion 
mass is clearly unphysical.We expect /x c ~ 
because the proton is the lightest state with non- 
zero baryon number. From this we are drawn to 
two possible conclusions. The first possibility is 
that the quenched approximation is at fault so 
that it is strictly necessary to consider the com- 
plex action of full QCD at finite chemical poten- 
tial. Secondly, there could be intrinsic problems 
in the lattice formulation of fermions (possibly as- 
sociated with fermion doubling) and chemical po- 
tential which would survive an unquenched treat- 
ment. 



Further studies |p[ of the quenched theory on 
larger lattices found similar behaviour. How- 
ever, a recent study |QJ|] of the quenched the- 
ory which measured the condensate and the pion 
and nucleon masses, did find that, at intermedi- 
ate and strong coupling, the theory is sensitive 
to the baryon mass but that it is pathological for 
M > ^- ' 

From an analytical study of the eigenvalues of 
the fermionic propagator matrix GibbsQ| con- 
cluded that the eigenmodes of the propagator 
matrix (defined below) relate to the mass spec- 
trum of the theory. The Gibbs argument, re- 
quires the calculation of the hadronic spectrum 
on replicated lattices, i.e. lattices, with periodic 
boundary conditions on the gauge fields, which 
have been strung together d times in the time di- 
rection. He considered the limit d — ► oo in order 
to replace finite sums with contour integrations. 
The procedure is justifiable at zero temperature. 
The expression obtained for the inverse of the du- 
plicated fermion matrix G(t\,t2) is 

G{t u t 2 ) = Y J A a \ t J- t * (1) 

where A a are the amplitudes (which can be re- 
lated to the eigenvectors of the propagator ma- 
trix) and A a are the associated eigenvalues. The 
form of the inverse shows that the exponential de- 
cay at large separation is controlled by the eigen- 
values A a : thus we see that the eigenvalue spec- 
trum calculated on isolated configurations should 
contain poles in correspondence to the physical 
masses. In particular, the smallest mass state can 
clearly be extracted from the lowest eigenvalues. 
This state is obtained by the squaring the propa- 
gator matrix, and defines the pion mass in QCD. 
This implies that in the quenched theory < ipip > 
is controlled by the pion mass. 

It has also been suggested by several authors 
that the coincidence of the onset of the chiral 
symmetry restoration with one half of the pion 
mass might only have been a numerical accident, 
the correct relationship being ^ onS et = ln ^ L — A 
where A is the contribution of the nuclear binding 
energy. If this scenario were true then the prob- 
lems with finite density simulations would not 
be too serious. This conjecture has been tested 
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pjj. However this work confirmed the onset at 

Recent work by Stephanov || using the random 
matrix method has provided some insight into 
this particular problem by demonstrating that at 
nonzero y, the quenched theory is not a simple 
rif — > limit of a theory with n / quarks but it is 
the nf — > limit of a theory with n / quarks and 
n./ conjugate quarks and therefore inappropriate 
to QCD.This explanation implies that the phase 
of the determinant is important for a simulation 
of full QCD at finite [i. In the quenched model 
the early onset for the number density would also 
correspond to the restoration of chiral symmetry 
because of the simultaneous occurrence of quarks 
and conjugate quarks in the system whereas the 
inclusion of dynamical fermions would result in 
a rearrangement of the eigenvalues such that the 
chiral symmetry would be restored at fi c ~ 
as expected. 

3. Methods for Simulating QCD at fi ^ 
with Dynamical Fermions 

The pathologies of quenched studies of the chi- 
ral phase transition at finite density indicate that 
a correct implementation of dynamical quarks in 
the simulations is essential to the physics of the 
model. 

3.1. The Glasgow Method 

The problem caused by the non-Hermitian na- 
ture of the fermion matrix at fj, ^ can be cir- 
cumvented by a method which involves the ex- 
pansion of the grand-canonical partition function 
(GCPF) in powers of the fugacity variable e M / T . 

The GCPF can be expressed as an ensemble 
average of the fermionic determinant normalised 
with respect to the fermionic determinant at /i = 
0. 



Z = 



J[dU}[dUi}\M(ri\ ( 



Sg[U,U*] 



f[dU][dW]\M(iJi = 0)\e- s ^ u < ut 1 



(2) 



where S g is the pure gauge action and M is the 
Kogut-Susskind fermion matrix. 

This means that the probability measure is 
well-defined and standard Hybrid Monte-Carlo 



algorithms can be applied. The GCPF is then 
given by the ensemble average 



/ |m(m)I \ 

\|MOu = 0)|A o 



(3) 



The Kogut-Susskind fermion matrix has the 
structure: 



2%M = G + 2im + Ve^ + F f e^. 



(4) 



in which G contains all the spacelike gauge links 
and the bare quark mass and V (V') all the for- 
ward (backward) timelike links. 

It then follows that, on a lattice of size n^nt, 
the determinant of M is given 



(5) 



\2iM\ = e 3fin > nt \P - e~i*\ 

where the propagator matrix P is given by 



P 



-{G + 2im)V V 



-V 







with its inverse given by 



P- 1 = 



-ft 

-V\G + 2im) 



(6) 



(7) 



From the structure of the propagator matrix 
and its inverse, it is easily shown that its eigenval- 
ues have the symmetry that if A is an eigenvalue 
then so is . In addition, since the matrix V 
factors, the eigenvalues have a Z(nt) symmetry. 
As a consequence of this symmetry, the charac- 
teristic polynomial for P is a polynomial in e pn * 
with 6n 3 + 1 complex coefficients. Since the prob- 
ability for a given configuration of gauge fields 
to appear in our ensemble is equal to that of its 
complex conjugate their imaginary parts will av- 
erage to zero. We impose this symmetry in the 
simulations described below. As a consequence, 
the characteristic polynomial from a configura- 
tion (averaged with its complex conjugate) is in- 
variant under fj, — ► — [i and has the form: 



3n/ 

E 

n— — 3n s 3 



(8) 
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The chemical potential dependence of the 
GCPF is then given explicitly by 

3n s 3 



-3n s 3 



3ri ^ 



E 



-0„-n M )/T 



(9) 



with 



: < b\ n \ > 



By measuring the e„ in the simulations we can 
readily obtain the fermion number density, 
p =< n > and the energy density, E =< e„ >. 

The relative values of the canonical partition 
functions can be used to characterize the proper- 
ties of the system. In particular the relative value 
of the triality non-zero coefficients give an indica- 
tion of whether the system is in the confined or 
deconfined regime. In the confined sector the en- 
semble average of the triality non-zero coefficients 
should tend to zero as the statistics increase. We 
see strong evidence for this behaviour is seen in 
the lattice QCD simulations described below. 

We also consider a description of the system 
in terms of the canonical partition functions for 
fixed particle number. This technique allows us to 
include only those triality zero coefficients which, 
at the end of the simulation, are positive and have 
a reasonable error. The chemical potential as a 
function of the baryon number density is obtained 
from the local derivative of the energy, e„ , of the 
state with n fermions with respect to n. 

"M " Sf£ (10) 

where p is the fermion density, g^j. 

3.2. Lee- Yang Zeros 

According to the theorems of Lee and Yang||] 
the phase structure of a simulated system can be 
determined by studying the zeros of the GCPF. 
The zeros correspond to singularities in the ther- 
modynamic potential. We study the zeros of the 
GCPF in the complex chemical potential (or fu- 
gacity plane) which are given by Z(p) = 0. Ac- 
cording to the Lee- Yang theorem as the system 
volume tends to infinity, the zeros will converge 
towards the critical value in the physical domain 
(real axis) of the chemical potential. 



3.3. Alternative Methods 

The Glasgow method involves an ensemble gen- 
erated at fj, = 0. This may introduce systematic 
errors due to insufficient overlap when measuring 
observables at ji significantly greater than the on- 
set of non-zero number density. An alternative 
method involves generating the ensemble either 
with respect to the absolute value of \M(u) or 
with respect to the modulus of its real part Jl0|, lT| . 
In these simulations it is necessary to measure an 
additional observable related to the average sign 
of the determinant. For example, updating with 
respect to the modulus gives the observable O as 



O 



| AT | 



(e 



(11) 



\M\ 



where (f> is the phase of the determinant. Simula- 
tions at a fixed quark mass using these methods 
are well behaved at low and high chemical poten- 
tials but for ^ < fi <~ 2£ the phase <f> fluctu- 
ates violently from configuration to configuration 
and < e 1 ^ > becomes very small and is not mea- 
surable. This is the region of p in which critical 
behaviour is expected. 

Gocksh|l2| used a spectral density method by 
binning the phase <fi and measuring p(<p). On a 2 4 
lattice, his results for p c are in broad agreement 
with meanfield. 

4. The Static Quark Model 



Blum, Hetrick and Toussaint 13 have recently 
studied numerical simulations of lattice QCD in 
the limit that the quark mass and chemical po- 
tential are simultaneously made large. In that 
limit, the quark mass and chemical potential ap- 
pear only in the ratio (2ma/e Aia ) n * and the prop- 
agator matrix becomes 



P = 



-2imV V 



-V 







(12) 



The corresponding fermion determinant is com- 
plex but trivial to evaluate which allows genera- 
tion of very high statistics in their measurements 
and determine < e 1 ^ > to sufficient accuracy. 

For quenched QCD the high temperature tran- 
sition is first order and they expected this be- 
haviour to extend into the interior of the T — 
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/i phase diagram. However their simulations 
showed that this transition becomes a smooth 
crossover at very small density (possibly for any 
nonzero density) and that, at low enough tem- 
perature, chiral symmetry remains broken at all 
densities. 

Of course, as the authors point out, it is not at 
all clear that the static approximation has any- 
thing to do with real QCD. However, it is rele- 
vant that unexpected results follow in this simple 
model. 

5. Lattice QCD at Strong Coupling 

The strong coupling limit of QCD provides 
us with an important testing ground for lattice 
Monte-Carlo simulations at finite density because 
we can compare our results with the analytic pre- 
dictions of mean-field theory. No such analytic 
predictions are available in the scaling region (i.e. 
intermediate coupling). At infinite coupling the 
mean-field method predicts a first order transi- 
tion at chemical potential fj, m f where 



Mm/ = - sinh 1 (X r) 



V 



(d- l)r 



where 



Ao 



rV2 



1 + (d- 1)V - 1 



1/2 



(13) 



(14) 



for a lattice with d space-time dimensions where 
r = n t /n s . 

The mean-field baryon and pion masses are 
given by 



TOfl = hi 



1 3 
-C J 

2 



1 r 

-c 6 



l + -(c 2 -2d) 



In 



+ \/(c 2 -2d) + -(c 2 -2d) 2 



where 



c = (m + y/2d + to 2 ) 



(15) 

(16) 
(17) 

(18) 



Note that [x m f does not correspond to the mean- 
field baryon threshold tob/3. 

It has been argued that the discrepancy is 
related to the binding energy of nuclear mat- 
ter which is large at infinite coupling. Bilic et 
al.jl4| demonstrated explicitly that 1/g 2 correc- 
tions diminish the discrepancy between /z m / and 
the mass of the lightest baryonic state divided 
by N (the number of colours). The mean-field 
pion thresholds and \x m f are not well separated 
for bare quark masses ma > 0.5. 

Monte-Carlo SU(3) simulations [ft5| using the 
monomcr-dimer algorithm whereby the partition 
function is represented in terms of monomers, 
dimers and baryonic loops estimate a critical 
chemical potential, fi c which is in agreement with 
fjL m f- This algorithm involves integrating out the 
gauge fields exactly and subsequently integrating 
the quark fields. Hybrid Monte-Carlo techniques 
integrate the quark fields prior to the gauge fields. 
The monomer-dimer algorithm cannot be used to 
explore small bare quark masses to < 0.1 and is 
restricted to the infinite coupling regime where 
the quarks are point-like. 

5.1. Strong Coupling using the Glasgow 
Method 

Consider again the expression for the GCPF 



3n s 



E 



-(e n -nn)/T 



(19) 



n— — 3n s 3 



The zeros of this polynomial are the Lee- Yang ze- 
ros of Z in the complex fugacity plane (or equiv- 
alently the complex e M plane). In the dynamical 
theory these zeros play a role analogous to that of 
the eigenvalues of P in the quenched theory when 
calculating the number density. 
Z can be rewritten as 



Z 



,3n>/T 



n ( z ~ 



(20) 



i— l,6n s 



where z = e~ M / T and the a« are the zeros of the 
averaged characteristic polynomial for P. 
On each configuration 



det (M) 



_ p 3«j M /T 



n ( z -^ 



(21) 



=1.6n s 
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where the A^ are the eigenvalues of P. Since 
we can interpret the zeros 



detM O=0) 



H=0 



of the partition function as the proper ensemble 
average of the eigenvalues of the fermionic propa- 
gator matrix. One obvious constraint on the dis- 
tribution of the zeros is that in the confined sector 
we expect to see a Z(3) symmetry arising from the 
triality non-zero coefficients averaging to zero so 
that only canonical partition functions containing 
multiples of three quarks contribute to the ther- 
modynamics. The Z(3) symmetry constrains the 
phases of the zeros but not their modulus. 

We calculate the fermion number density, n q 
from 



T d\n{Z) 
V dfj, 



(22) 



Hence in the quenched theory the number den- 
sity, n q , is given by 



1+ < 



1 



3F 

V 



3V ^-f \i~ z 

-3V 



> . 



(23) 



where the A^ are the eigenvalues of P nt . Using 
the generalized boundary conditions arising from 
considering a configuration on a replicated lattice 
oneQ can convert the sum to a contour integral 
and obtain, for each configuration 



n d = 1/V h 



(24) 



1< | A; \<e^l T 



In the dynamical theory the number density, 
rid, is given by 



nd 



1 



1 

W 



3V 

E 

-3V 



(25) 



where the cti are the zeros of the averaged charac- 
teristic polynomial for P. Performing the contour 
integation as above gives 



n d = 1/V 



E 

K|ai|<ef/ T 



1. 



(26) 



Thus if we consider plotting the zeros of the 
GCPF in the complex e^-plane, any discontinu- 
ities in the number density will be associated with 
circular bands of increased density of zeros. 



5.2. Results at Strong Coupling 

We first consider, in the light of the above dis- 
cussion concerning the roles of the eigenvalues 
and the zeros of the GCPF, their distributions 
in the complex fi and planes. 




Figure 1. Histograms of the real parts of the 
eigenvalues and of the real part of the zeros in 
the complex /i-plane at strong coupling and bare 
quark mass, ma — 0.08 on a 6 4 lattice. 



Fig contains a histogram of the distribution 
of the positive real parts of the eigenvalues in 
the complex fi plane. The distribution is ob- 
tained from 450 configurations on a 6 4 lattice at 
a bare quark mass of 0.08. There is a clear sig- 
nal that the number density becomes non-zero 
at \i ~ 0.3 which is consistent with m n /2. This 
is the quenched onset \i. The lattice is filled at 
H = 0.95. 
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The analogous distribution of the zeros of the 
GCPF is also shown in Fig [3. These zeros are 
those of the polynomial obtained by averaging the 
characteristic polynomials of the propagator ma- 
trix over the same 450 configurations. There is 
still a clear signal that the number density be- 
comes non-zero at the same onset fi as that of 
the quenched theory but a strong signal has devel- 
oped via an intermediate peak in the distribution 
which is absent from that of the eigenvalues. As 
argued above, this band of increased density can 
be associated with a discontinuity in the number 
density. 

This difference in distribution between the 
eigenvalues and the zeros is found at all other 
bare quark masses. Fig. |^ and Fig. || show our 
results at quark mass m = 0.05 while Fig. | and 
Fig. H show results at m — 0.5. 

Figs. Hand ^show clearly the banded structure 
of the zeros distribution in the e M plane for two 
different lattice sizes and quark masses. 

Fig. |^ compares the distribution of the real 
part of the zeros with the number density (on a 
6 4 lattice with m — 0.1). Again the histogram 
for the zeros shows three distinct peaks: y, Q ~ 0.3 
corresponding to the onset of net non-zero quark 
density; /i c ~ 0.7 corresponding to the small 
discontinuity in the number density and the ex- 
pected critical chemical potential; \i s ~ 1.0 corre- 
sponding to the lattice saturation point. Compar- 
ison of the two plots in this figure shows that the 
derivative of the number density correlates well 
with the density of states (frequency histogram). 

Fig. ^ gives an overview of the location of 
fj, a and fi c for six quark masses in the range 
ma = 0.05 to 1.5. The scaling of the central peak 
of the histogram is in remarkable agreement with 
the mean-field prediction calculated form eqn.13 
and is consistent with the critical /i predicted by 
the monomer-dimer simulations. The peak corre- 
sponding to fi a is consistent with the mean-field 
pion threshold for the smaller quark masses but 
significantly lower for ma — 1.0, 1.5. 



6. Lattice QCD at Intermediate Coupling 



histogram of eigenvalues 

ma=0.05 (650 configs at strong coupling) 



10000.0 




Re (log WLT) 



Figure 2. Histogram of the real parts of the 
eigenvalues in the complex /i-plane at strong cou- 
pling and bare quark mass, ma — 0.05 on a 6 4 
lattice. 



6.1. Z(3) Tunnelling 

In the pure gauge theory we have N c (= 3) 
equivalent vacua related by Z(N C ) rotations. 
Tunnelling between the different Z(3) vacua is 
much more probable in the confined sector than 
in the deconfincd sector. Since the pure-gauge 
action as well as the integration measure are in- 
variant under the Z(3) transformation, the GCPF 
can also be written as: 

J[dU][dW]detM^ + z 3 M)e-^ f ] 
[P> J[dU}[dW}detM(ii = 0)e- s ^ U ' U ^ { ' 

Averaging over the three Z(3) vacua would elim- 
inate the triality non-zero coefficients and, since 
Z(0) = 1, the sum of the triality zero coefficients 
should tend to 1. As Fig. ^0| shows, at — 5.1 and 
m = 0.01, we obtained a clear signal for Z(3) tun- 
nelling. The sum of the triality zero coefficients 
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Histogram of distribution of zeros 

ma =0.05 



150.0 



0.0 



-2.0 -1.0 



0.0 
Re n 



1.0 2.0 



histogram of evs 

ma=0.5 (200 configs at strong coupling) 



4000.0 - 



2000.0 - 



0.0 
0.60 



0.80 0.90 
Re (logWLT) 



1.00 1.10 



Figure 3. Histogram of the real parts of the 
zeros in the complex /i-plane at strong coupling 
and bare quark mass, ma — 0.05 on a 6 4 lattice. 



Figure 4. Histogram of the real parts of the 
eigenvalues in the complex /z-plane at strong cou- 
pling and bare quark mass, ma — 0.5 on a 6 4 
lattice. 



shows a strong tendency to average to 1 whilst 
the sum of the triality one (and two) coefficients 
tends to zero. This is consistent with simulating 
in the confined phase at /i = as expected at this 
coupling and quark mass. 

6.2. Simulations on 6 4 and 8 4 lattices at in- 
termediate coupling 

We have investigated [^6| the chemical poten- 
tial (the onset /i) required to make the level with 
3 quarks equally probable with that of the zero 
quark level, i.e. fionset = (^3 — eo)/3- This ad-hoc 
definition takes into account that it is these en- 
ergy levels which are most accurately determined 
and allows errors to be estimated directly. 

Fig. [ll] shows the dependence of the onset 
H on the quark mass and on the square root of 
the quark mass on a 6 lattice at gauge coupling 



(3 — 5.1. As at strong coupling, there is a strong 
signal that this onset fj, is dependent on the mass 
of a Goldstone boson for quark mass m > 0.01. 
(On a 6 4 lattice, at this /3, the system becomes 
'deconfined' at some m < 0.01. We have checked 
this behaviour by determining its Lee- Yang ze- 
ros in the complex m-plane and they are complex 
with non-zero real part < 0.01.) 

Fig. [12] shows the distribution of the zeros 
found on a 8 4 lattice at the same (3 and quark 
mass 0.01 in the e M plane. The zeros associated 
with the onset /z lie on the arc of a circle (this cir- 
cular pattern is also observed in the strong cou- 
pling results described above). Again the distri- 
bution is consistent with a number density be- 
coming non-zero at a small /i. Also we see a cir- 
cular band of zeros developing at 0.7 < /1 < 0.9. 
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Histogram of zeros distribution 



ma=0.5 



300.0 



200.0 



0.0 
Rejt 



E 

9- - 



Re(exp(mu)) 



Figure 6. Zeros in the plane for m = .095 , 
6 4 lattice . The critical line is the thin line inside 
the denser region e M = e Mc 



Figure 5. Histogram of the real parts of the 
zeros in the complex /i-plane at strong coupling 
and bare quark mass, ma — 0.5 on a 6 4 lattice. 



We believe that this is a signal for the physical 
Me- 



7. The U(l) Gross-Neveu Model 

7.1. The Lattice Formulation 

In the above simulations of QCD, the initial 
onset of the critical behaviour at non-zero \i was 
controlled by the pion mass rather than that of 
the lightest baryon. To check if this is an artifact 
arising from measuring our observables on an en- 
semble generated at [i = we have investigated 
the above methods in a study of the Gross-Neveu 
model at finite density. The lattice formulation 
of this model is described in ref. (L7| 

The lattice action for the bosonized Gross- 



Neveu model with U(l) chiral symmetry is 

i=l *-x,y 

X 

+ is(x) tt(x) 

<x,x> 



(28) 



Here, Xi an d Xi are complex Grassmann- valued 
staggered fermion fields defined on the lattice 
sites, the auxiliary scalar and pseudoscalar fields 
a and n arc defined on the dual lattice sites, and 
the symbol < x, x > denotes the set of 8 dual 
sites x adjacent to the direct lattice site x. Nf is 
the number of physical fermion species. 

The fermion kinetic operator M is given by 



M 



x.y 



J y,x-0 



10 




Re(exp(mu)) 



Figure 7. Zeros in the e M plane for to = .1 
lattice. 



1 



+ 

2 

i/=l,2 

+ TnSy iX , 



J y,x-\-v 



Oy,x — 



(29) 



where to is the bare fermion mass, /x is the chem- 
ical potential, and r) v {x) are the Kawamoto-Smit 
phases (— l) a;rH The symbol denotes 

the alternating phase (— i^o+xi+m . 

This lattice model (at non-zero lattice spac- 
ing) has the symmetry: U(iV//4)y ®U(./V//4)y ® 
U(1) J 4- It is the U(1)a symmetry which is bro- 
ken, either spontaneously by the dynamics of the 
system, or explicitly by a bare fermion mass. We 
simulated the Nf = 12 model corresponding to 
three lattice species. 

The expansion of this GCPF is very similar to 
that described above for standard QCD. At finite 
density the Dirac fermion matrices M and M are 
given by: 

2iM xy {^) = Y xy + G xy + V xy e» + V^e"" 
-2iM xy (f,) = Y2 y + G xy + V X ye» + V2 v e-». 



where Y = 2i{m q + \Y, <XtX> {cr{x) + ie-K{x)))8 X y. 
The determinants of these fermion matrices are 
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Figure 8. The number density at quark mass 
m q = .1 (upper) on a 6 4 lattice at strong cou- 
pling. The behaviour of the number density is 
reflected in the histogram of the real parts of the 
zeros in the complex /i-plane shown in the lower 
portion of the figure. 



related to that of the propagator matrix 
-GV-YV V 



P = 



-V 







(30) 



by det(2iM) = e 3 ^™' det(P - e~») and 
det(2iM) = e 3 ^"* det((p- 1 ) t - e~»). 

As for standard QCD, determination of the 
eigenvalues of P nt gives the expansion for the 
GCPF: 



2n a 



Z = J2 < Vl > e 
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Figure 9. Summary of the results for the critical 
point fi c , and current onset fi . [i c follows the pre- 
diction of the mean field analysis of ref. [4] . The 
onset is close to half the pion mass at small mass, 
and lower for m q > .5 . The data for the pion 
mass are from ref. [7] 



2n B 2 

= £ e -(«»-»M)/T (31) 

n=—2n B 2 

7.2. Results 

This model has a massless pion and its fermion 
determinant is real and non-negative even when 
/i / 0. Therefore standard Monte-Carlo algo- 
rithms can be used to study its critical behaviour 
as a function of \i. 

In particular, at an inverse four-fermion cou- 
pling of 0.5 and bare fermion mass m = 0.01, 
the theory has a critical chemical potential /i c = 
0.725(25) at which the chiral symmetry is re- 
stored. Note, the Gross-Neveu model is non- 
confining and the gauge fields do not enter into 
the dynamics. 

We performed simulations developing the en- 
semble at fi = to confirm that Eqn. 5 does 
give an extrapolation which predicts the correct 
critical behaviour. The number density for a 16 3 
lattice at the above coupling and fermion mass 
is shown in Fig. [l3|. The chemical potential at 
which the number density, < n >, becomes non - 
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Figure 10. The HMC time evolution of the ex- 
pansion coefficients on an 8 4 lattice at f3 = 5.1 
and ma = 0.01. The behaviour of the triality 
2 coefficients is similar to that of the triality 1 
coefficients. 

zero, Ho, is correctly predicted. Comparison with 
the exact simulation \\n\ (in which generation of 
the ensemble was at the fi at which the measure- 
ment was made) shows that at \x > fi Q the num- 
ber density is underestimated and the disconti- 
nuity in < n > is absent. The Lee- Yang zeros in 
the complex fi plane are the zeros of Eqn. 5 and 
should signal the critical fi associated with the 
chiral transition transition. They are shown in 
Fig. m There is a line formed by 6 zeros, two 
of which are isolated from the rest of the dis- 
tribution. The line intersects the real fi axis at 
/i = 0.72, in agreement with the critical /i found 
by the exact method. 

Hence we find that the onset \x found via the 
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Figure 11. The quark mass dependence of the 
onset fj, ((3 = 5.1 on a 6 4 lattice). 

number density is not controlled by the Goldstone 
boson, and the critical fi is given by the distribu- 
tion of the Lee- Yang zeros. 

Fig. 13 also shows the number density obtained 
from an ensemble generated at an update fj, = 0.8, 
greater than [i c . In the chiral limit it is consis- 
tent with that of a free gas of massless fermions 
whereas the number density from the [i = up- 
date is that of a free gas of massive fermions. The 
exact results show the transition between these 
two phases at /i c . 

8. Conclusions 

Simulation of standard lattice QCD at finite 
density seems to be plagued by the problems of 
the naive quenched version of the theory, namely 



Figure 12. The Lee- Yang zeros on a 8 4 lattice in 
the complex e** plane at = 5.1 and bare quark 
mass 0.01. Only one segment is shown but the 
intrinsic Z(n t — 8) of Figs. 4 and 5. also applies 
here. 



the chemical potential at which the fcrmion num- 
ber begins to differ from zero is controlled by 
a Goldstone pion rather than by the lightest 
baryon. 

Clearly questions must be answered: 
Why do we continue to observe an onset chem- 
ical potential which scales with the pion mass? Is 
the band of zeros at higher fi a signal of the tran- 
sition in the number density? Is this band asso- 
ciated with the deconfinement and chiral transi- 
tions? Do these two transitions occur at the same 

The success of the monomer dimer approach 
(of use only at strong coupling) may be due to 
analytic integration first over the gauge fields, 
and then the fermion fields, giving an effec- 
tive action with a reduced sign problem. Stan- 
dard Hybrid Monte-Carlo algorithms integrate 
the fermion fields prior to the gauge fields. We 
are now attempting to understand why at strong 
coupling, we find an onset chemical potential, fi , 
in addition to the expected signal at fi c . 

Incomplete cancellation could arise from either 
limited statistics or from the reweighting proce- 
dure itself. Using an ensemble created at /x = 0, 
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Figure 13. The fcrmion number density vs. the 
chemical potential, ^t, on a 16 3 lattice at an in- 
verse four-fermion coupling 0.5 and fermion mass 
0.01 



Figure 14. The Lee- Yang zeros in the complex 
\i plane for the Gross-Neveu model at an inverse 
four-fermion coupling 0.5 and fermion mass 0.01 
on a 16 3 lattice 



together with reweighting with respect to the 
fermion determinant at \i = may not give suffi- 
cient overlap with an ensemble generated at fi ^ 
to correctly describe the true physics there. Al- 
though a limited increase in statistics revealed 
that the peak in the histogram of zeros corre- 
sponding to [i Q grew proportionally with the peak 
at \x c it is still possible that the signal for ji 
will be cancelled out in the limit of infinite statis- 
tics once the Z(3) invariance has been completely 
achieved. 

Another solution may be to approach QCD by 
adding a perturbatively irrelevant four-fermion 
interaction term to the action and to study the 
critical nature of that theory as this additional 
term is weakened. This approach permits simu- 
lations directly in the chiral limit. 
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